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Abstract In this paper, the velocity field and the time dependent shear stress of Maxwell fluid with 
Caputo fractional derivatives in an infinite long circular cylinder of radius R is discussed. The 
motion in the fluid is produced by the circular cylinder. The fluid is initially at rest and at time 
t = 0 + , cylinder begins to oscillate with the velocity fsimot , due to time dependent shear stress acting 
on the cylinder tangentially. The hybrid technique used in this paper for the solution of the problem 
has less computational efforts and time cost as compared to other commonly used methods. The 
obtained solutions are in transformed domain, which are expressed in terms of modified Bessel 
functions I 0 (j and /](•). The inverse Laplace transformation has been calculated numerically by 
using MATLAB package. The semi analytical solutions for Maxwell fluid with fractional deriva¬ 
tives are reduced to the similar solutions for Newtonian and ordinary Maxwell fluids as limiting 
cases. In the end, numerical simulations have been performed to analyze the behavior of fractional 
parameter a, kinematic viscosity v, relaxation time X, radius of the circular cylinder R and dynamic 
viscosity /i on our obtained solutions of velocity field and shear stress. 

© 2017 Faculty of Engineering, Alexandria University. Production and hosting by Elsevier B.V. This is an 
open access article under the CC BY-NC-ND license (http://creativecommons.Org/licenses/by-nc-nd/4.0/). 


1. Introduction 

Many viscoelastic fluids have applications in industry and 
engineering. It contains the extrusion of polymer fluids, sus¬ 
pension and colloidal solutions, biomagnetic fluids, food 
stuffs, exotic lubricants and cooling of the metallic plate. 
Newtonian and Navier-Stokes view is based on constitutive 
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equations of the fluids. Non-Newtonian fluids are the fluids 
that cannot be modeled by Navier-Stokes equations. The 
non-Newtonian fluid describes all those fluids, due to which 
shear rate is not directly proportional to shear stress. 

The book of Chanderashker [1] is considered standard for 
the exact solutions of non-Newtonian fluids in cylindrical 
domains. Makris and Constantinou [2] discovered that it is 
not possible for Maxwell fluid to get sufficient experimental 
data due to range of its different frequencies. Friedrich [3] 
found that there is a relationship between governing equations 
involving fractional derivatives and molecular theories. Makris 
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et al. [4] presented that when the ordinary Maxwell fluid has 
been replaced by Maxwell fluid with fractional calculus, a good 
fit of experimental data is obtained. The constitutive equations 
of non-Newtonian and rate type fluids are studied by Rajago- 
pal [5 and differential type fluids have been studied by Dunn 
and Rajagopal [6], Unsteady pressure flow of a classical Max¬ 
well fluid in a cylinder is discussed by Rahaman and Ramkis- 
soon [7], In particular, the predictions of a fractional Maxwell 
model with the linear viscoelastic date in glass transition and a 
relaxation zeros is given by Palade et al. 8]. Wood [9] studied 
the helical movement of an Oldroyd-B fluid due to the com¬ 
bined action of a constant pressure gradient and rotating cylin¬ 
ders. Zhu et al. [10] studied the characteristics of the shear 
stress and the velocity field of an ordinary Maxwell fluid. Some 
exact solutions with respect to different flows of non- 
Newtonian fluids with fractional derivatives in cylindrical 
domain have been presented in [11-14]. 

Peristaltic flow of fractional Maxwell fluid model through a 
medium examined by Tripathi et al. 15]. Peristaltic transport 
of fractional Maxwell fluid model in uniform tubes was studied 
by Tripathi [16]. The motion of a fractional Maxwell fluid 
inside two walls perpendicular to a plate is discussed by Vieru 
et al. [17]. Yang and Zhu [18] studied start up motion of a vis¬ 
coelastic fluid in a pipe with fractional calculus approach. The 
movement of a viscoelastic fluid between two side perpendicu¬ 
lar walls with vibrating shear stress is discussed by Fetecau 
et al. [19], Oscillating flows of Newtonian and Maxwell fluid 
with slip and no slip effects can be found in the references 
[20-22]. Tripathi [23] discussed peristaltic transport of a Max¬ 
well model in a channel where Tripathi and Anwar [24] found 
numerical solution of oscillating peristaltic flow of viscoelastic 
fluids of a generalized Maxwell model through a porous med¬ 
ium. Tripathi and Anwar [25] investigated peristaltic flow of 
viscoelastic fluid with a slip condition. 

Exact analytical solutions for a longitudinal flow of a frac¬ 
tional Maxwell fluid between two coaxial cylinders is discussed 
by Awan et al. [26]. The unsteady flow of a Maxwell fluid with 
fractional derivatives in a circular cylinder moving with a non¬ 
linear velocity is discussed by Athar et al. [27]. Soret and 
Dufour effects on unsteady flow of an Oldroyd-B fluid with 
boundary conditions discussed by Ashraf et al. [28]. Some 
interesting and recent results for the solution of viscoelastic flu¬ 
ids with fractional derivatives have been presented in [29-33]. 

In several industrial, engineering and biological transport 
processes, oscillating motion of viscoelastic fluids in a tube 
or channel is a fundamental flow dominion confront. The 
examples of oscillating flow in biological systems are the 
quasi-periodic blood flow in the cardiovascular system, urody- 
namic transport in the human ureter and movement of food 
bolus in the gastrointestinal tract. The examples of oscillating 
flows in industrial systems include waste conveyance systems 
employing roller pumps and finger pumps. As such fractional 
models have been studied extensively in recent years in 
biomedical transport problems. 

This paper deals with the semi analytical solutions for the 
unsteady flow of a Maxwell fluid with fractional derivatives. 
Initially the fluid in the circular cylinder of radius R is at rest 
and at time t = 0 + , because of the shear stress, it starts to oscil¬ 
late longitudinally with the velocity fsinwt. The semi analytical 
solutions are procured with the tool of Laplace transformation 


and modified Bessel equation. The inverse Laplace transforma¬ 
tion has been achieved through MATLAB. The obtained solu¬ 
tions are in series form in terms of modified Bessel functions 
/o(-) and / 1 (-) and meet all initial and boundary conditions. 
We obtained the known solutions corresponding to a Newto¬ 
nian fluid for 2, a —> 0 and ordinary Maxwell fluid for a —> 1 
performing the same motion, as the limiting cases of our gen¬ 
eral solutions. According to the best of our knowledge, no 
work has been done to solve oscillating Maxwell fluid models 
with this technique. A table is included in the paper, which 
shows the good agreement of our obtained results of inverse 
Laplace transform, with two other existing numerical algo¬ 
rithms, named as Stehfest’s [41] and Tzou’s [42], Finally, some 
graphs are presented to show the impact of physical 
parameter. 

It is worthy pointing out that the integral transform meth¬ 
ods previously used for finding the exact solutions of such kind 
of models have some flaws. For example the graphical solu¬ 
tions obtained by Hankel transformation in [34] does not sat¬ 
isfy its initial and boundary conditions. Recently, Abdullah 
et al. [35] used the numerical Laplace method in an efficient 
way to solve fractional Maxwell model. Sheng et al. [36] found 
that the numerical inverse Laplace algorithms are reliable for 
fractional differential equations. Tahir et al. [37], Raza et al. 
[38], Tong et al. [39] and Jiang et al. [40] used the numerical 
inverse Laplace algorithm successfully to solve fluid models. 

2. Basic governing equations 


The velocity field G and the extra-stress T for the movement of 
fluid [9] are considered as 

G = G (r,t)=g(r,t)e z , T = T (r,t), (1) 

where e- is the unit vector of the cylindrical coordinate system 
(r, 9 , r) along the r-direction. Moreover, when the fluid starts 
to move, we have 

G(r, 0) = 0, T(r, 0) = 0. (2) 


The governing equations, due to such motions for Maxwell 
fluids [9] are 

d 


f ) = 9 


d \ dg(r, t) 

x+x jt)^r 


dg(r, t) 
dr ’ 


(3) 

—■ l ~ 

+ 

<0 | co 

> 

| 9g(r, t) 

(4) 

1 dr ’ 


where ((r, t) = T r0 (r , t) is the only nontrivial shear stress, /r is 
the dynamics viscosity, v is the kinematics viscosity and X is 
relaxation time of the motion. The Caputo fractional differen¬ 
tial operator [9] is defined as 


D e ,g(t) 


1 £ r> g(0 J r 

r(i- e ) dt Jo (<-{)« 

<feU 

dt ’ 


0 ^ q < 1 ; 
6 = 1, 


(5) 


where T(-) denotes the Gamma function. By using the Caputo 
fractional differential operator Z)“ instead of inner time deriva¬ 
tives in Eqs. (3) and (4), the equations related to the Maxwell 
fluid having fractional derivatives can be obtained as 


(1 + FD\%r, t)=n 


dg{r, t) 
dr 


( 6 ) 
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(l + /TD“) 


dg{r, t ) 

dt 


, a 1 

- V[ dr + 7- 


9g(r, t) 
dr 


(7) 


3. Oscillatory flow through a circular cylinder 


Consider a Maxwell fluid with fractional derivatives which is 
initially at rest in an infinitely long circular cylinder of radius 
R. In the presence of shear stress, at/= 0 + , the cylinder begins 
to oscillate due to shear stress with the velocity fsinmt. The 
fluid is gently stirred due to the shear with the velocity of the 
form (1),. 

The initial and boundary conditions of the flow are 


g(r,0) = 


dg(r, t) 


dt 


= 0, ((»', 0) = 0; re[0,iJ], 


(1 + ^)C(M)U = /i : 


dg{r, t) 


dr 


= fsinmt , 


( 8 ) 


(9) 


where/is a constant and co is the angular frequency. The Eqs. 
(6) and (7) involve fractional derivatives, these equations are 
solved by using the tool of Laplace transformation and modi¬ 
fied Bessel equation. 


4. Calculation of the velocity field 


Taking the Laplace transformation of Eqs. (7) and (9), we have 


9 2 g(r,s) 1 dg(r,s) 

dr 2 r dr 


(v + 2V +1 ) 

v 


g{r,s) = 0 , 


( 10 ) 


dg(r, s) = fco 
dr r=R li(s 2 + co 2 ) ’ 


( 11 ) 


where s is the Laplace transform parameter. Eqs. (10) and (11) 
can be written as 


<?g{r,s) 1 dg(r,s) 

dr 2 r dr 


a(s)g(r, s) = 0, 


( 12 ) 


0 b(s) 

1 y/aisjliiRy/aisj)’ 

substituting the value of C i in Eq. (16), we get 
y/a{sjli(Ry/a{sj) 


(17) 


(18) 


Eq. (18) is in a complicated form of modified Bessel func¬ 
tions of first and second kind. For the solution of Eq. (18), it 
is very difficult to find the inverse Laplace transform tradition¬ 
ally. To solve this problem, we used some numerical packages 
to obtain the inverse Laplace transform. Here, we have found 
the inverse Laplace transform numerically through MATLAB. 
Our technique is referred as semi analytical because the analyt¬ 
ical and numerical methods are used in this paper. 


5. Calculation of the shear stress 


Taking Laplace transform of Eq. (6), we have 


C(r,s) 


__f__ dg{r, s) 

(1+2V) dr 


Differentiating Eq. (18) with respect to r, we have 

dg{r,s) _ h( s 7 i(''\/«W) 
dr ’ IfR^ffs)) 


(19) 


( 20 ) 


Substituting the above Eq. (20) into Eq. (19), we obtain 


C(r,s) = b(s) 


g h{rs/a(s)) 

(1+AV) IfR^dfs))' 


( 21 ) 


After taking inverse Laplace transform of Eq. (21), we 
obtain the solution of shear stress. The expression of Eq. 
(21) is of the same type of Eq. (18 ). Here, again we have found 
the inverse Laplace transform numerically through MATLAB. 


6. Limiting cases 


dg{r,s) 

dr 

where 


= b{s), 


,21 ' 


(13) 


(s + X i“ +1 ) fco 

a(s) = - -and b(s) = —f — 

' v w ji(s 2 + ( 

By introducing variable transformation m = r-Ja(s) in Eq. 
(12), we get 

m 2 + m^- — (in 2 — 0 2 )g = 0. (14) 

dm 2 dm ' V 


The general solution of this modified Bessel equation is 
given by 

g(m,s) = CJoim) + C 2 K 0 (m ), (15) 


where C\ and C 2 are constants, Io(m) and K 0 (m) are the mod¬ 
ified Bessel functions of the first and second kind respectively. 
In order to have a finite solution at m = 0 (r = 0), C 2 must be 
zero. Then Eq. (15) becomes 

g(m,s) = C t I 0 (m). (16) 


Using the boundary condition given by Eq. (13) with Eq. 
(16), we obtain 


We can take the following limiting cases of our general 
solutions. 

6.1. Ordinary maxwell fluid 

Using a —> 1 in Eqs. (18) and (21), we acquire the transformed 
solutions of the velocity field g 0 M( r , t) and the shear stress 
CoAc( r ,t) respectively equivalent to ordinary Maxwell fluid. 

6.2. Newtonian fluid 

When we take 0 in Eqs. (18) and (21), the transformed 

solutions of the velocity field g N (r, t) and the shear stress 
Cj v (r, t) equivalent to Newtonian fluid are procured 
respectively. 

7. Results and discussion 

In order to provide the validation of results, we presented in 
Table 1, numerical results of fluid velocity, obtained with 
MATLAB program and with other two numerical algorithms, 
namely the Stehfest’s algorithm [41] and Tzou’s algorithm [42]. 
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Table 1 Comparison between different numerical algorithms. 


r g(r, t ) (MATLAB) 

[43] 

g(r, t) (Stehfest’s) 
[41] 

g(r, t ) (Tzou’s) 
[42] 

0 

0 

0 

0 

0.05 

0.449567 

0.446226 

0.447221 

0.1 

0.451195 

0.447853 

0.448848 

0.15 

0.453910 

0.450566 

0.451561 

0.2 

0.457715 

0.454368 

0.455364 

0.25 

0.462615 

0.459264 

0.460261 

0.3 

0.468615 

0.465259 

0.466257 

0.35 

0.475722 

0.472361 

0.473360 

0.4 

0.483944 

0.480577 

0.481577 

0.45 

0.493291 

0.489916 

0.490918 

0.5 

0.503772 

0.500389 

0.501392 



According to Stehfest’s algorithm, the inverse Laplace trans¬ 
form is given by 


, In 2A 

,0 = — 

1 k= 1 


Ini 


r, 


V k = 





(f -mu 


ftm _ 

-\)\(k-j)\(2j-k)V 


( 22 ) 


where N is the number of the expansion terms and must be an 
even number (N = 16 leads to a very good precision). 

Tzou’s algorithm is based on the Riemann-sum approxima¬ 
tion. In this method the inverse Laplace is given by 



(23) 

where Re(-) is the real part, i is the imaginary unit and N\ is a 
natural number. The values of inverse Laplace transform 
obtained with MATLAB, Eqs. (22) and (23) are given in the 
following table. 

It is observed that the results with different numerical algo¬ 
rithms has a good agreement among them. 

Finally, we have plotted some graphs for the velocity field 
and the shear stress of the fluid by using Eqs. (18) and (21) 
respectively, to see the effect of various material parameters 
on our results. Fig. 2(a) and (b) depicts r with shear stress 
f(r, /) and velocity field g(r , t ) for various values of t at 
H = 15, /= 10, A = 5, R = 0.1, co = 5, 
a. = 0.5, v = 0.035754 and it is observed that the shear stress 
C(r, t) and velocity field g(r, t) is an increasing function. Other 
graphs have been plotted against the values of t. Fig. 3(a) and 
(b) depicts t with shear stress and velocity field g(r, t) 

for various values of relaxation time A = 5 at 
H= 15, /= 10, co = 5, R = 0.1, a = 0.5, v = 0.035754. We 
conclude that the shear stress and the velocity field are decreas- 


7 



°0 0.02 0.04 0.06 0.08 0.1 


r 

(a) Shear stress 



(b) Velocity 


Fig. 2 Shear stress £(r, t) and the velocity field g(r, t) graphs of the fractional Maxwell fluid by using material parameters 
/i = 15, /= 10, A = 5, R = 0.1, co = 5, a = 0.5, v = 0.035754 and various values of t. 
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Fig. 3 Shear stress £(r, t) and the velocity field g(r, t) graphs of the fractional Maxwell fluid by using material parameters 
g = 15, /= 10, co = 5, R = 0.1, a = 0.5, v = 0.035754 and various values of 2. 



Fig. 4 Shear stress f(r, t) and the velocity field g(r, t) graphs of the fractional Maxwell fluid by using material parameters 
g = 15, /= 10, 2 = 5, co = 5, 7? = 0.1, a = 0.5 and various values of v. 


ing function to relaxation time 2. Fig. 4(a) and (b) depicts t with 
shear stress f(r, t) and velocity field g(r, t) for various values of 
kinematic viscosity v at // = 15, / = 10, 2=5, co = 5, 
R = 0.1, a = 0.5 which show that the kinematics viscosity v is 
increasing function to both the shear stress and the velocity 
field. Fig. 5(a) and 5(b) depicts t with shear stress £(r, t) and 
velocity field g(r, t) for various values of the radius of circular 
cylinder R at /c = 15, /= 10, 2 = 5, co = 5, v = 0.035754. 
a. — 0.5 which show that shear stress and velocity field are 
decreasing function to the radius of circular cylinder R. Fig. 6 
(a) and (b) depicts t with shear stress £(r, t) and velocity field 
g(r, t) for various values of the fractional parameter a at 
g = 15, /= 10, 2=5, R = 0.1, ca= 5, v = 0.035754 which 


show that shear stress and velocity field are decreasing function 
to the fractional parameter a. Fig. 7(a) and (b) depicts t with 
shear stress f(r, t) and velocity field g(r. t) for various values 
of the angular velocity co at /c= 15, /= 10, 2 = 5, 
v = 0.035754, a = 0.5, R = 0.1 which show that the shear 
stress and velocity field are increasing function to angular 
velocity co. Fig. 8(a) depict t with velocity field g(r. t) for various 
values of the dynamics viscosity v at n = 15, /= 10, 2 = 5, 
co = 5, a = 0.5, = 0.1 which shows that dynamics viscosity 

is decreasing function to velocity field. Fig. 8(b) depict the com¬ 
parison of velocity field g(r, t) between Maxwell fluid with frac¬ 
tional derivatives, ordinary Maxwell and Newtonian fluid at 
H= 15, /= 10, 2 = 5, 
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(a) Shear stress ( b ) Velocity 


Fig. 5 Shear stress £(r, t) and the velocity field g(r, t) graphs of the fractional Maxwell fluid by using material parameters 
/i = 15, / =10, 2 = 5, oj = 5, v = 0.035754, a = 0.5 and various values of R. 



time time 


(a) Shear stress (b) Velocity 

Fig. 6 Shear stress £(r, t) and the velocity field g(r, t) graphs of the fractional Maxwell fluid by using material parameters 
/i = 15, /= 10, 2=5, R = 0.1, co = 5, v = 0.035754 and various values of a. 


v = 0.035754, co = 5, a = 0.5, R = 0.1. It is seen from Fig. 8 
(b) that Newtonian fluid is swiftest and Maxwell fluid with frac¬ 
tional derivatives is the slowest. 

8. Conclusions 

In this paper, the flow of Maxwell fluid governed by the differ¬ 
ential equations with Caputo fractional derivative are studied. 
The flow domain is the inner of a circular cylinder and flow is 
generated by the longitudinal stress-force given on the cylinder 
surface. By applying the Laplace transform with respect to the 
time variable ?, the semi analytical solutions for the fluid veloc¬ 
ity and longitudinal shear stress are obtained in terms of the 
modified Bessel functions of first kind 7 0 (-) and /((•). Since 


the Laplace transforms of the velocity and of the shear stress 
are enough complicated, we have obtained the inverse Laplace 
transforms by means of the numerical procedures. Firstly, we 
used a MATLAB numerical code. In order to provide valida¬ 
tion of results, we have used other two numerical algorithms, 
namely the Stehfest’s algorithm and Tzou’s algorithm. 

As shown in Table 1, we found a good agreement between 
the results obtained for inverse Laplace transform with three 
numerical methods. It is important to observe that the fluid 
layers situated close cylinder surface have a significant motion, 
while the fluid situated in the central area of the cylinder has a 
very slow motion. The fluids modeled with fractional deriva¬ 
tives flow faster than the ordinary fluid. When the fractional 
parameter decreases, the fluid velocity increases. Shear stress 
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time 

(a) Shear stress 


time 

(b) Velocity 


Fig. 7 Shear stress f(r, t) and the velocity field g(r, t) graphs of the fractional Maxwell fluid by using material parameters 

/y = 15, y = 10, 2=5, R = 0.1, a = 0.5, v = 0.035754 and various values of to. 




time 

(a) Velocity 


time 

(b) Shear Stress 


Fig. 8 Velocity field g(r, t) graph with various values of n and comparison graph of velocities of fractional Maxwell, ordinary Maxwell 
and Newtonian fluid by using material parameters /r = 15, /= 10, 2 = 5, v = 0.035754, to = 5, a = 0.5, R = 0.1. 


has the behavior similar with velocity; therefore, it is increasing 
when the fractional parameter decreases. Physical parameters 
t. co and v are directly proportional to both velocity field 
and shear stress and parameters R, a, /t and 2 are opposite 
to both velocity field and shear stress. 
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